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Massive planets in FU Orionis discs: implications for 
thermal instability models 
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ABSTRACT 

FU Orionis are young stellar objects undergoing episodes of enhanced luminosity, 
which are generally ascribed to a sudden increase of mass accretion rate in the sur- 
rounding protostellar disc. Models invoking a thermal instability in the disc are able 
to reproduce many features of the outburst, but cannot explain the rapid rise time- 
scale observed in ma ny cases. Here we explore the possibility (originally suggested by 
iClarke fc Sverlll996j) that the thermal instability is triggered away from the disc inner 
edge (at a distance of « IORq from the central protostar) due to the presence of a 
massive planet embedded in the disc. We have constructed simple, one-dimensional, 
time-dependent models of the disc evolution, taking into account both the interaction 
between the disc and the planet, and the thermal evolution of the disc. We are indeed 
able to reproduce rapid rise outbursts (with rise time-scale ~ 1 yr) , with a planet mass 
M s = 10 — 15Mj up itor- We show that the luminosity and the duration of the outbursts 
are increasing functions of planet mass. We also show that the inward migration of the 
planet is significantly slowed once it reaches the radius where it is able to trigger FU 
Orionis outbursts, thus suggesting that a single planet may be involved in triggering 
several outbursts. 

Key words: accretion, accretion discs - instabilities - stars: formation - stars: pre- 
main-sequcnce - stars: variables (FU Orionis) 



1 INTRODUCTION 

FU Orionis objects are a small class of Young Stellar Ob- 
jects undergoing violent and probably recurrent outbursts, 
during which they can increase their bolometric luminosity 
by tw o to three orders of magnitude JHartmann fc KenvorJ 
1996). These outbursts are generally attributed to a sud- 
den increase of the mass accretion rate (which can reach 
10 _4 MQ/yr) in the disc of an otherwise "normal" T Tauri 
star. This interpretation is suggested by a number of ob - 
servations: (i) in one case (i.e. V1057 Cyg, iHerbia ll977T) 
a pre-outburst stellar spectrum is available, and it shows 
typical features of a T Tauri star; (n) the Spectral Energy 
Distribution (SED) after the o utburst is well described in 
terms of accretion dis c SEDs IIKenvon fc H artmann 1991; 
lLodato fc Bertir]|200lf) (Hi) optical and near infrared line 
profiles are usually double-peaked, as expect ed if they are 
produced in a differentially rotating d isc l|Keny on etal 
I198SI: lLodato fc BertirJl200a see however Eerbig et alJl2005 
for a dissenting view). 

The light curves of the three best studied FU Orionis 
objects (i.e., FU Ori itself, V1515 Cyg, and V1057 Cyg) 
show remarkable differences between each other. The rise 
timescale of FU Ori and of V1057 Cyg is very short (of the 



order of lyr), while that of V1515 Cyg is definitely longer 
(trise ~ 20 yrs). On the other hand, while FU Ori and V1515 
Cyg have a very long decay timescale (idccay ~ 50 — 100 yrs), 
V1057 Cyg decays much faster (td Gca y ~ 10 yrs). 

Several different mechanisms have been suggested to 
trigger the outbur st: these include a tidal interaction with 
a companion star JBonnell fc BastieiJll992l). a gravitationa l 
instability in the outer, massive disc dArmitagee^ay|2001), 
and a viscous-thermal instability (IBell fc Lira 11994 . here- 
after BL94) in a disc fed at a high enough mass accre- 
tion rate from a surrounding envelope. In particular, the 
latter model (which had been extensively studied in or- 
der t o explain the dwarf novae outbu rsts in binary sys- 
tems, iMever fc Mever-Hofmeisterlll98lf) has been discussed 
in some details and is considered to be the most promising 
model for the outburst. 

Despite its ability to reproduce ma ny feat ures of the 
outburst, the thermal instability model feL94l) has found 
some difficulties in explaining the details of the FU Orionis 
light curves. In fact, a thermal instability is triggered when 
the disc surface density becomes larger than a critical value 
Ea, dependent on the relevant opaciti es in the quiescent 
phase (see Sect. [5] below). In the IBL94 models this occurs 
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first very close to the inner edge of the disc, so that the 
instability propag ates inside-out, producing a slowly rising 
luminosity. |BL2J indeed find rise timescales > 10 yrs, which 
may be compatible with the light curve of V1515 Cyg, but 
are too long to explain the rapid rise of FU Ori and of V1057 
Cyg. On the other hand, a rapid rise can be achieved if the 
instability is first triggered at a large radius, s o tha t the 
instability propa gates outside-in. IClarke et al.l 1119901) and 
iBell et all il995T> have constructed models where the insta- 
bility is triggered by introducing some ad hoc density pertur- 
bations at R > 1O7?0 from the central object, and succeeded 
in obtaining rapid rise outbursts. 

An alternative (and self-co nsistent) way to obtai n rapid 
rise outbursts was proposed bv lClarke fc Sv er (1996). They 
showed that the presence of a companion with mass larger 
than the local disc mass (for example, a massive planet) em- 
bedded in the disc can substantially modify the disc struc- 
ture, banking up the disc material just outside the location 
of the companion. In this way, the banked up disc surface 
density may eventually become larger than the critical value 
and the instability can be triggered at the radius correspond- 
ing to the position of the companion. 

IClarke fc Armitagel J2003I) have recently suggested that 
a planet embedded in the disc of an FU Orionis object would 
lead to a clear spectroscopic signature in the form of a peri- 
odic modulation of the double-peaked line profiles observed 
in these systems, with periods corresponding to the orbital 
frequency of the planet. Interestingly, periodic modulations 
in the line profiles have indeed been observed in a long term 
monitoring campaign o f the rapid rise FU Orionis objects 
V1057 Cyg and FU Ori jHerbig et all2003l) . The persistence 
of the feature (with a period of w 3 days) in successive ob- 
serving series argues against its origin as a transient disc 
surface feature, and is consistent with the hypothesis of an 
embedded hot Jupiter. 

In this paper we develop a simple, one-dimensional, 
time-dependen t outb urst model following the suggestion of 
IClarke fc SveJ 1119961) that the outburst is triggered by the 
presence of a planet embedded in the disc. We take into ac- 
count the details of the planet-disc interaction and of the 
thermal evolution of the disc. In this way we are able to 
follow at the same time both the orbital evolution of the 
planet, the banking up of disc material outside the planet 
position and the details of the outburst phase. We indeed 
obtain rapid rise outbursts, with the rise timescale and the 
location of the radius at which the disc is first triggered into 
outburst depending mainly on the planet mass. [Note that in 
this paper we are interested in the dynamical effects of the 
presence of a companion with small mass ratio embedded 
in the disc, to which we will refer as a "planet" , even if we 
will allow its mass to be above the nominal minimum mass 
for deuterium burning, sometimes used to define planets as 
opposed to brown dwarfs.] 

The paper is organised as follows: in Section[5]we intro- 
duce the thermal instability model and discuss the specific 
disc model that we have adopted here. In Section [3] we dis- 
cuss the details of the planet-disc interaction. In Section 2] 
we show the results of our outburst simulations. In Section 
|^| we discuss our results and in Section |S| we draw our con- 
clusions. 



2 THERMAL INSTABILITY MODEL 

2.1 General features 

The time evolution of an axisymmetric, thin, Keplerian ac- 
cretion disc (in which, for the moment, we do not include 
any effects coming from the tidal interaction with an embed- 
ded p lanet) is determined by the diffusion equation JPringlel 

1981): 



9E 3d ( 1/2 d 1/2 
m = RdR{ R 8R {R A 



(1) 



where R is the cylindrical radius, E is the disc surface den- 
sity, and n — vYj = M /2,tt is the mass flux at radius R, 
where v is the viscosity coefficient. The viscosity is usually 
described in terms of some empirical prescri ption, such as 
the a prescription (Sha kura fc Sunvaevlll973h . according to 
which v — ctc 3 H, where c s is the disc thermal speed and H is 
the disc thickness. The value of \i can be obtained from the 
solution of the energy balance equation for the disc, once a 
prescription for the opacity (which determines the relevant 
radiative cooling term) is given. In this paper, we will use 
the following sim plified form of the t ime-dependent energy 
balance equation l|Pringle et al.l ll986'l: 



dfj, 



t-thcrm 



/i9E 
E dt' 



(2) 



where £ t hcrm = a _1 S7 _1 is the thermal time-scale, and Q 
is the angular velocity of the disc. In equation J2J, /j, eq is 
the value of /j, at thermal equilibrium, and is determined 
by the specific choice for the opacity. In equation @ we 
have neglected some terms, namely the radial radiative flux 
(which is going to be of the order of H/R with respect to 
the other terms, and therefore negligible for a thin disc) and 
the advective terms (i.e. terms in the form u • V, where u is 
the fluid velocity). The simplified description adopted here 
is able to ca pture the main prop erties of the outburst (see, 
for example. iPringle et al.lfl986l) . even if some quantitative 
result may actually be influenced by the terms that we have 
neglected (see Section |S] below). 

A given equilibrium solution is unstable if <9/i oq /<9E < 0. 
This happens, for example, when the temperature of the 
disc is close to 10 4 K and hydrogen becomes partially ionised 
(BL94). The behaviour of a thermally unstable disc is best 
understood in the E-/i plane at a fixed radius, where the 
equilibrium curve can have a characteristic S-shape, in which 
two stable solutions are connected by an intermediate un- 
stable branch (see Fig. 0. The two critical values for the 
density, Ea and Eb , represent the maximum disc density on 
the lower branch and the minimum disc density on the up- 
per branch, respectively. If, for some reason, the disc density 
becomes greater than Ea the disc cannot stay in the lower 
stable branch and jumps (on a thermal timescale) to the up- 
per branch, characterised by a much larger accretion rate. 
As a result of the enhanced accretion rate, the surface den- 
sity will decrease and eventually become smaller than Eb, 
and the disc will go back to the lower quiescent branch. Of 
course, the process described above refers to an isolated disc 
annulus, but when the instability is triggered somewhere in 
the disc it pro pagates to the a djacent annuli with a typical 
velocity » aCs JLin et al.lll985l) . As a result, when the insta- 
bility propagates inside-out, since the thermal velocity of the 
disc decreases with increasing radius, it will slow down dur- 
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Figure 1. Typical S-curve in the S-/x plane at a given radius. 
Sa is the maximum disc density on the lower branch, Eb is the 
minimum disc density on the upper branch. 



ing the propagation, producing a slowly rising light curve. 
In addition, viscous diffusion always results in a higher mass 
fraction evolving inward and this, together with the fact that 
in general Ea decreases at small radii (see below) makes the 
inward front propagation faster. 

The time spent by the disc on the upper branch can 
be estimated simply as the viscous time-scale at the outer 
radius where the instability propagates: 



t v 



R\m 



>CRlim) ' 



(3) 



2.2 Reference disc model and numerical method 

As we have seen in the previous Section, to set up a thermal 
instability model we need to specify the details of the equi- 
librium curves (the S-shaped curve in Fig. at each radius. 
In this work we will use approximate paramctcrisations of 
the detailed vertical structure calculations by BL94: 
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where Ei ow and Ehigh are the cold and hot stable branches 
of the equilibrium curve, respectively, and where Ea and 
Eb are the critical values of the density (see Fig.Q. We will 



also assume that the unstable branch is a straight line in the 
plane logE-log/i, connecting the points (Ea, ma) and (Eb, 
^b) (see Fig. 0. 

In thermal in stability models of dwarf novae outbursts 
(see lLasotall200ll for a review) it has become customary to 
adopt a larger value of the viscosity parameter a in the upper 
branch with respect to the lower branch, in order to correctly 
reproduce the light curves of the outburst. This might have 
some physical justification in the scenario in which the main 
mechanism for angular momentum transport in the disc is 
MHD turbulence: in fact, in the lower branch the disc has a 
lower degree of ionisation with respect to the higher branch, 
so that magnetic instabilities are expected to be less effective 
in promoting accretion. In this work, we will conform to this 
practice and will adopt two different values for a in the two 
equilibrium branches. Our reference values are cvi ow = 10 -4 , 
Qhigh = 10~ 3 , but we have also explored different choices 
for a (see below). The transition from the low to the high 
value of a is followed smoothly when the disc moves from 
the lower to the upper branch and viceversa. 

We have solved the coupled equations JTJ and (2) using 
an explicit scheme, with a radial grid equally spaced in the 
variable R 1 ' 2 . The inner and outer disc radii are taken to be 
Rin = 3Rq and R ou t = 150Rq, respectively. We have taken 
a free flow inner boundary condition and a reflecting outer 
boundary condition. 

In all our outburst models we have fed the disc at an 
accretion rate M[ n kept constant with time. This is simply 
implemented by adding a source term in equation Q. This 
term is taken to be vanishing outside a narrow radial range 
close to the outer boundary of the disc. 

To save computational time, if all the disc lies on the 
lower branch (and if E is everywhere smaller than Ea), 
rather than solving equation (J5J, we simply assume that \x is 
equal to its equilibrium value on the lower branch, which is 
equivalent to the requirement that the density at each point 
satisfies equation J1J . In this way the limiting factor for the 
time-steps in our code is simply the stability condition for a 
viscous flow, which generally gives a time-step much longer 
than the thermal time-scale on which we should evolve the 
code if we also solved equation @ . 



2.3 Simulated inside-out outbursts 

As a test of our numerical scheme, we have first simulated 
so me insi de-out outbursts of the same kind of those obtained 
by|BL9J, without a ny perturb ation (neither artificially im- 
posed, as in Clarke et al. 1990, nor self-consistently obtained 
through the presence of a planet). 

Fig. [5] shows the light curves obtained in this way, for 
an imposed acc retion rate Mi n — 3 10 -6 MQ/yr (so as to be 
consistent with BL94|). The central object mass is assumed 
to be Mi, = Mq (we will use this assumption throughout 
the paper). We are indeed able to obtain repetitive out- 
bursts, with a duration of ~ 100 yr a nd a re currence time 
of ~ 2000 yr, in fair agreement with IBL94I . The peak of 
the bolometric luminosity is I/ pca k ~ 10 2 Lq, and the mass 
accretion rate during the outburst is M ~ 6 10 Mq/jt. 
The rise in luminosity is characterised by an initially rapid 
increase, which eventually slows down so that the maximum 
is achieved after ~ 10 yr, in agreement with the expecta- 
tions from an inside-out outburst. The outburst is indeed 
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an inside-out outburst, the radius where the instability is 
first triggered being ~ 4Rq. The outer radius to which the 
instability propagates in this ca se was Ry lm ~ 25Rq, again 
in agreement with the results of[BL94. 

In other simulations, not shown here, we have also been 
able to reproduce the variation of the outburst features with 
respect to parameters variations (dif ferent c hoices of a, dif- 
ferent values of M 1B ) as described in IBL94 . The agreement 
that we obtain is remarkable, considering that our outburst 
m odel is significantly simplified with respect to the one used 
bv |BL94| . In particular, (i) while w e use simple parameter- 
isations for the equilibrium curves, IBL94I solve explicitly in 
much greater detail the vertical disc structure and the radia- 
tive transfe r, inclu ding also non-thermal-equilibrium compu- 
tations; (m)|BL9J retain all the advective terms and the ra- 
dial energy fluxes in their thermal equation, while we neglect 
them. Our results show that the main qualitative features 
of the outburst are not significantly affected by these more 
refined calculations. 



3 PLANET-DISC INTERACTION 

In this Section we will leave for a while the subject of ther- 
mal instability in discs and will concentrate on the second 
basic ingredient of our model, i.e. the interaction between 
an accretion disc and an embedded satellite. It is worth 
noting that the calculations described in this Section can 
be applied also to similar systems on a much larger scale, 
such as the case of a binary of two supermassive black holes 
jNataraian fe Armitagell2002J) . 



3.1 Generalities: gap formation, planet migration 

The topic of disc-satellite interaction and the con- 
nected issue of satellite mig ration has been explor ed 
in a number of papers (e .g., iLin fc Papaloizoul Il979al lbi 
iGoldreich fc Tremainelll98 0i. It is well known that a mas- 
sive enough satellite is going to open up an annular gap in 
the disc at its orbital radius as a result of the gravitational 
torques exerted between the satellite and the disc. This an- 
gular momentum exchange will also cause the satellite to 
migrate to smaller radii, in the so-called Type II migration 
(as opposed to the much faster Type I migration which oc- 
curs when the planet mass in not high enough to open up a 
gap in the disc). 

In this paper we will use a simplified expression for the 
tidal torque, obtained in the "impulse approximation", by 
analogy with the process of dynamical friction. In this case, 
equation Q should be modified as follows: 



dt RdR\ 8R ( N ) RdR\Q 
where the specific tidal torque A is given by: 



A = q 2 n 2 R 2 ( ^ 
p 



A = -q 2 Q 2 R 2 



R 



R> R* 



R< Rs 



(8) 



(9) 



M B /Mi, is the mass ratio between the satellite and the cen- 
tral object, and p = R — R s . This simplified form of the spe- 
cific torque is th e same as commonly used by many investi- 
gator s (see, e.g.. lArmitage et alJl2001i lArmitaee fc BonnelJ 
2002). This formalism, although simplified, has been shown 
to provide results i n agreement with m ore sophisticated nu- 
merical treatments (ITrilling et alll QQSI. In addition, for the 
massive planets that we are considering here, the form of 
the torque plays a minor role in the resulting dynamics (at 
least during the quiescent phase of the limit cycle), since 
the satellite is massive enough to produce a clear gap. In 
such Type II migration, the evolution of the planet is in- 
sensitive of the form of the torque because the disc surface 
density profile self-adjusts so as to lock the planet orbital 
migration to the evolution of the edge of the gap. The spe- 
cific form of the torque will be more important during the 
outburst phase, when the gap is going to be filled. In this 
case, estimates of the torques (and of th e timescale for Type 
I migration) based on the linear theory (ITanaka e t al. 2002) 
might not be applicable for a number of reasons: (i) the in- 
teraction with the high mass planets that we consider here 
might be non-linear; (n) these estimates usually neglect the 
torques coming from inside the Roche radius of the plan- 



ets, which might be important BcU^^^iL[[2003); (ra) nm- 



iii) 
yvE> 



away migration (as described bv lMasset fc Papalo izou 2003) 
might influence significantly the results, when the disc mass 
is high. Most 2D and 3D si mulations of planet-disc interac- 
tion i n this embedded case <|Lubow et alJll999l : iBate et alJ 
2003: Mas set fc Papaloizoull2003i) are not useful to our pur- 
poses, since they generally explore the case where the disc 
is a much colder, T Tauri-like disc, in which case a planet 
that does not open up a gap needs to have a small mass 
(less than I Mj up jt G r). Detailed numerical simulations of high 
mass planets embedded in hot discs would be needed to 
clearly assess this problem (incidentally, such simulations 
would also have the chance of better resolving the region 
within the Roche lobe of the planet, since this is going to be 
relatively large). Such analysis is, however, beyond the im- 
mediate goals of the present paper. In the absence of clear 
results, we have therefore preferred to use the above formu- 
lation for the torque also during the outburst. 

For numerical reasons, we have smoothed the torque 
term when R ~ _R S , where the torque would have a singular- 
ity (see equation JUJ). Wc have used the same smoothing pre- 
scription as in lSver fc Cla rke (1993) and lLin fc PapaloizovJ 
(1986), so that the torque is smoothed when \R — R B \ < 
max[_ff, Rh], where H is the disc thickness and Rn is the 
size of the Hill sphere of the planet (note that, for the disc 
models and planetary masses considered here, during quies- 
cence generally Rn > H, while during the outburst H and 
Ru are comparable and of the order of w 0.I-R). 

The back reaction of the disc on the planet orbital mo- 
tion can be obtained by simply imposing angular momentum 
conservation, that can be written in the form: 



-(M s n s R 2 s ) = 



2nRAZdR, 



(10) 



where the integral is taken over the whole disc surface. 

The behaviour of the coupled disc-satellite system de- 
pends on two dimensionless parameters: 



In equation ©, R s is the radial position of the satellite, q 



. MR 2 2 
A= q 2 , 



(11) 
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Figure 2. Light curves of simple inside-out outbursts. In this case, the imposed mass accretion rate was M- ln = 3 10 Mq/yi. The right 
panel shows a blow up of the outburst light curve. 



measures the relative strength of the second and first terms 
on the right-hand-side of equation J5J, while: 



B = 



4tyT,R J 



(12) 



gives a measure of the magnitude of the right-hand-side of 
equation 1101 . 

The parameter A h as a simple physical inter pretation. 
In fact, it can be shown llLin fc Papaloizoui ri979b) that A = 
(A/_R) 3 , where A is the gap width. In order to open up a 
gap, the gravitational effect of the satellite has to overcome 
the pressure of the disc and th e viscosity, that both act as 
to counteract the gap opening. Takcuchi ct al. ( 1996) have 
shown that a gap is going to be open when q > 2a 1 ' 2 (H/R) 2 . 

The second important parameter, B, represents the rel- 
ative local disc mass with respect to the satellite. If the 
satellite is more massive than the disc (i.e. B <C 1) its in- 
ertia will cause a rather slow migration, whereas in the op- 
posite case (B ~^> 1) the planet will migrate on the local 
viscous timescale, behaving like a representative fluid ele- 
ment of the disc. In the limit where B <C 1 the satellite will 
act like a dam for the flow, preventing accretion through the 

J ap, and will cause th e upstream surface density to bank up 
Sver fc Clarke! 1995). In general, B is an increasing func- 
tion of radius. We can therefore define a radius Rb , at which 
B = 1, so that for R ^$> Rb the planet will be light compared 
to the disc and it will simply migrate with the accretion flow. 
For example, for the quiescent surface density profile that 
we adopt here (equation @), we have: 



Rb « 0.62i? G 



io- 



AL 



M 



J y-L Jupiter 
0.8 



lO- 6 M /yr 



(13) 



3.2 Simulated planet migration 

ISver fc Clarke! lll995D and llvanov et all <1999h have found 
analytical solutions for the orbital evolution of the planet 
and for the banking up of surface density in the disc in 
the approximation in which the gap width is assumed to 
be very small and equation JHJ is substituted by an appro- 
priate boundary condition at the planet position. In these 
computations the inner disc (i.e. the disc inside the satellite 
orbit) is assumed to be empty, and initially the outer disc 
is assumed to have the surface density corresponding to an 
equilibrium disc accreting at a constant rate. 

As we have already mentioned, the evolution of the sys- 
tem (concerning the satellite-disc evolution) is "scale-free" 
and is determined only by the two parameters A and B, 
and by the functional relation between E and /i, which is 
normally taken to be in the form E oc fi a R b (cf. equations 
(0J and ©). If initially the disc is in a steady state at some 
given accretion rate M, the satellite evolution, according to 
llvanov et alJ (119991) . is given by: 



R s = R s0 {l->yB (t/t ) d Y 



(14) 



where i? s o is the initial position of the planet, and 7 and d 
are two dimensionless parameter of order of unity, depen- 
dent only on a and b. Bo is the initial value of B (in the 
unperturbed disc) and to is the viscous timescale in the un- 
perturbed disc at R = Rsq: 



t = 



4fi B 2 o 
Su(R B o) ' 



(15) 



In this work we have decided to tak e a different ap- 
proach with respect to llvanov et alJ Bl999T> and to solve nu- 
merically equations JHJ and 11U1 . In fact, (i) in the outburst 
simulations that we will present in Section 2] we will not 
have an initially steady state disc, but we will feed it from 
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Figure 3. Le ft pan el: Evolution of the satellite orbit in our numerical simulation (solid line) and according to the analytical solution of 
llvanov et alj Jl999I) (dotted line). Right panel: Initial surface density profile (solid line), compared with the surface density profile of the 
disc at three different times, i.e. £/£q = 1.7, 2.9, 4.1, respectively. 



outside at a given rate, and (u) one important goal of this 
work will also be to check what surface density profile will 
be left over in the inner disc after the outburst, and whether 
the satellite will still be able to open a gap at that stage. 
Therefore, we will need to treat in full details the density 
evolution of the disc at all radii. 

We have used the same explicit code discussed above 
in Section |^1 and the gravitationa l torque has been com- 
puted using the Lelevier method JPotteil 11973') . As men- 
tioned above, we have smoothed the torque term in equation 
lO . This smoothing will not affect the overall evolution of 
the satellite, but is going to prevent the gap from being fully 
open, so that a small leakage (with a magnitude of a few per- 
cent of the unperturbed accretion rate) from the outer to the 
inner disc is present in our simulations. Three-dimensional 
simulations of p lanet-disc interactions llLubow et alJll999t 
iBate et ah 2003) often show that a leakage through the gap 
is actually occurring, so that our model is not going to be un- 
realistic. This mass flow through the gap can be significant 
for low mass planets, but in the case we are considering here, 
with a high mass planet embedded in a cold, low- viscosity 
disc, a clear gap is going to be open, and the m ass flow 
through the gap is very small. IBrvden et al.l <l2000r) estimate 
the amount of mass flow through the gap as a function of 
the relevant disc parameters, finding that: 



dM 



■ 5 10~V 



I) 



5.2au 



200g/sec 

1/2 



M /yr. 



(16) 



During quiescence, the typical values for the disc density and 
thickness can be evaluated from equation Q. We find that 
at a typical radius of R « IORq , H/R sa 0.035. According to 



equation 116t . the mass flow through the gap is then going 
to be only « 2.5 10~ 10 MQ/yr, consistent with the numerical 
"leakage" that we allow in our simulations. 

We have checked our planet-disc interaction numerical 
scheme simulating the evolution of a sate llite in an initially 
steady-state disc (as in llvanov et alJll999l) . with no external 
feeding. The results are shown in Fig. [3] In this simulation 
the initial values of A and B were Aq = 0.1 and Bo — 0.1, 
respectively, and the equilibrium curve is given by equation 
@. With this choice for the surface density profile, B grows 
linearly with radius (see equation 1131 ). so that Rb — 10-RsO- 
The left panel of Fig. [3] shows the orbital evolution of the 
satellite (solid line), compared to the analytical prediction 
according to equation 114H (dotted line). The simulated mi- 
gration agrees reasonably well with the approximate analyt- 
ical expectation, being only slightly faster. (Note, however, 
that llvanov et alJll999l also numerically simulate the evolu- 
tion of a satellite, within the infinitesimal boundary layer 
model, and find as well a slightly faster evolution with re- 
spect to the analytical formula). 

The right panel in Fig. [3] shows the initial surface den- 
sity of the disc (solid line) , compared with the surface den- 
sity at three different times during the planet evolution. 
We can clearly see the formation of the gap and the in- 
ward migration of the sate llite (cf. analogous results in 
iNataraian fc ArmitageJ[2002l in the context of supermassive 
black hole binaries). Most importantly for our study, this 
plot shows clearly how the presence of the satellite is going 
to enhance the upstream surface density of the disc with re- 
spect to its initial value, and how the inner disc is rapidly 
depleted. In the next Section we will show how, as a con- 
sequence of this "banking up" of material, the upstream 
surface density can become larger than the critical value 
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for thermal instability Ea, thus triggering the outburst at a 
relatively large radius. 

As a final comment, we note that, under particular cir- 
cumstances, the migration of the planet can be significantly 
modified by the effect of corotation resonances, that we have 
not include d here, and th at might lead to a rapid "runaway" 
migration (Ma sset fc Papaloiz ou 2003). If the surface den- 
sity profile of the disc is extremely shallow (with E falling 
off with radius less rapidly than FC 1 ' 2 ) this type of migra- 
tion might even be outward. However, for these effects to 
be important, the gap should not be completely empty, so 
that the corotation region would co ntain enough mass to 
provid e a significant torque. Indeed, iMasset fc Papaloizou! 
(2003) only find runaway migration to occur in the transi- 
tion region between Type I and Type II migration, when the 
gap is partially open. For planet masses larger than lMj up it er 
th ey recover the stan dard Type II migration (see Fig. 12-14 
in lMasset fc Papaloizou! I2003T) . We therefore conclude that 
during the quiescent phase corotation torques are not going 
to significantly alter our picture. 



4 PLANET INDUCED OUTSIDE-IN 
OUTBURSTS 

We have obtained planet-generated outbursts in FU Orio- 
nis objects by combining the planet-disc interaction code 
described in Section [3] with the thermal evolution code de- 
scribed in Section |5] In what follows we have always used 
the equilibrium curves defined by equations 0-J7J. All our 
outburst simulations started with a relatively "empty" disc, 
i.e. a steady state accretion disc on the lower equilibrium 
branch accreting at relatively low rate, Mo = lO _lo M0/yr. 
This initially very light disc is fed at a constant rate dur- 
ing the simulation, as described in Section but, contrary 
to that case, now the matter which is added in the outer 
disc cannot flow freely to the smallest radii because of the 
banking up effect due to planet. 

In the next Subsection we will describe the details of 
the outburst generated in our "reference case" , in which the 
planet mass is M s — 10Mj up it O r, the disc is fed at Mi n = 
3 10~ 6 MQ/yr, and where aa ow = 10~ 4 and ahigh = 10~ 3 . 
The initial orbital radius of the planet is 40Rq . The values 
of the input mass accretion rate and of a are chosen so as to 
easily compare with the results of BL94. In Subsection l4.2l we 
will discuss the effects of changing the relevant parameters. 



4.1 Details of the reference case 

Fig. 2] shows the various phases of the disc evolution and of 
the outburst. In each plot the surface density of the disc is 
plotted at three different times, along with the critical val- 
ues Ea(-R) and Eb(-R), displayed with solid lines. Initially, 
the planet opens up a gap in the disc (as one would ex- 
pe ct from 3D simulat i ons in the case o f a massive planet, 
cf. iBrvden et alJll999t iBate et alJ l2003) and the inner disc 
is rapidly drained into the central star, leaving an almost 
empty inner "hole" , while the planet migrates inward (panel 
a). This first phase occurs on the slow viscous timescale. The 
three snapshots in panel (a) correspond to t « 23000, 28500, 
and 33000 yrs, respectively. When the surface density of the 
disc reaches Ea, the instability is triggered: in the reference 



case this happens at R as IORq ; the region of the disc where 
the instability is triggered moves to the upper branch and 
becomes much hotter, thus allowing pressure forces in the 
disc to overcome the tidal effect of th e planet: the gap is 
close d (cf. 3D numerical simulations bv lClarke fc Armitagd 
2003) and the instability front rapidly propagates through 
the inner disc (panel b). The typical timescale of the evolu- 
tion in this phase is the faster thermal timescale. The three 
snapshots in panel b refer to t m 35100 yr, and are sep- 
arated by one year between each other. Once the inward 
propagation front has reached the inner disc, the instabil- 
ity continues to propagate outside much more slowly, until 
the maximum extension of the propagation is reached: in the 
reference case this happens at Ry im as 43Rq . The mass accre- 
tion rate through the inner disc is M out ss 2.2 10 _4 MQ/yr 
(panel c, snapshots separated by « 5 yrs). The instability 
then slowly retreats from outside-in: in this phase the disc 
is moving from the upper branch to the lower branch and 
is still too hot for the planet to re-open the gap (panel d, 
snapshots separated by « 6 yrs). At later times, the disc 
cools down to the lower branch and the planet is now able 
to open the gap again (panel e, snapshot separation (a 5 yrs). 
The inner disc is again drained rapidly onto he star, while 
the outer disc is banked up again, being ready for another 
outburst to take place (panel f). This phase now occurs on 
the slower viscous timescale (snapshot separation in panel f: 
« 1000 yrs). 

Fig. |S] shows the light curve of the planet-induced out- 
burst compared to the light curve of the inside-out outburst 
described in Section^ All the parameters in these two mod- 
els are identical, except for the presence of the planet. 

The first thing to notice is that the outside-in outburst 
triggered by the planet indeed shows a faster rise time-scale. 
The maximum luminosity is reached in « 2 yrs, whereas in 
the inside-out outburst, the luminosity initially rises quickly 
but then slows down, reaching the maximum in ~ 7 yrs. 

A second important feature is that the peak lumi- 
nosity is much higher in the planet-triggered case, being 
ipcak ~ 43OZ/0 (in close agreement with the observed peak 
luminosity for FU Ori), to be compared to a maximum lu- 
minosity I/p Ca k ~ 1001/0 for the inside-out outburst. This 
can be easily understood by considering the equilibrium so- 
lutions, equations @-lJ7|l. In fact, during the outburst, the 
inner disc accretes at an approximately constant rate, de- 
termined roughly by the mass accretion rate in the higher 
branch corresponding to the radius where the instability is 
first triggered, where the surface density is Ea (i.e. to the 
accretion rate denoted in Fig. 0with /jaa). By combining 
equations Jjj) and © we can thus obtain an estimate of the 
outburst mass accretion rate: 



Mo 



3.5 10" 



Qhigh 

Q^low 



16/15 



/ _Rtrig_ 

\WR G 



M Q /yr, (17) 



where -R tr i g is the radius at which the instability is triggered. 
We then see that M ou t is a strongly increasing function of 
the trigger radius, so that an outside-in outburst is bound 
to be more violent than an inside-out one. In this respect is 
very interesting to notice that, among the three best studied 
FU Orionis objects, the two that show a rapid rise outburst 
(i.e. FU Ori and V1057 Cyg) are indeed characterised by a 
larger peak luminosity with respect to V1515 Cyg, which is 
characterised by a slower rise time. 
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Figure 4. Evolution of the disc surface density profile during the outburst in the reference case. The two solid lines show the critical 
surface densities Ea(-R) and Eb(-R)- (panel a: t = 20000 — 33000 yrs) The surface density is banked up upstream of the planet position; 
(panel b: t = 35176 — 35178 yrs) the outburst is triggered at R Rj 10-Rq and the inward propagation front of the instability rapidly makes 
its way to the inner edge of the disc; (panel c: t = 35180 — 35190 yrs) the outer propagation front slowly reaches the outer radius at 
which the instability can propagate; (panel d: t = 35200 — 35212 yrs) the instability retreats; (panel e: t = 35250 — 35260 yrs) the planet 
is able to re-open a gap; (panel /: t = 36000 — 38000 yrs) the disc is ready for another outburst. 
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Figure 5. Light curves of the planet-triggered outburst (solid 
line), compared to the inside-out outburst (dotted line). We have 
used here a linear scale for the luminosity to better emphasise 
the differences between the two cases. The off-set of the temporal 
axis has been opportunely chosen for display convenience. 




; yt 

Figure 6. Orbital evolution of the planet. During the first 3 10 4 
yrs the disc builds up, until the first outburst is triggered. The 
orbital evolution of the planet is significantly slowed down by the 
outburst mechanism. 



IClarke fc Sverl 1119961) predicted that the long term evo- 
lution of the orbit of the planet will be significantly slowed 
down by the outburst mechanism. This is because the sur- 
face density left over after the outburst will be small, and 
will therefore not be able to push the planet effectively. Once 
the surface density has increased enough to start moving the 
planet inward, then a new outburst will start and the density 
will drop down again. Fig. HJ shows the long term evolution 
of the planet position. During the first w 3 10 4 yrs the disc 
builds up, until the first outburst is triggered. Subsequently 
the evolutio n of the planet i s sign ificantly slowed down, as 
predicted by IClarke fc Sv cr (1996). Between t fs 3 10 4 yrs 
and t ~ 5 10 4 yrs the system undergoes ps 8 outbursts and 
the planet moves from R sa WRq to R w 3Rq . 

The slow inward migration also affects the light curves. 
In fact, as the planet moves inward, the instability is sub- 
sequently triggered at smaller radii. This causes the peak 
luminosity to become progressively smaller and the recur- 
rence time to decrease: initially, two consecutive outbursts 
are separated by w 4000 yrs, while by the end of the simu- 
lation the recurrence time has decreased to ~ 2000 yr (see 
Fig.0. 



;;; 
o 



(t/yr) 



Figure 7. Long term evolution of the luminosity in the "reference 
case". The peak luminosity and the recurrence time decrease for 
late outbursts. 



4.2 Effect of varying parameters 

4-2.1 Planet mass 

Table Q illustrates how the properties of the outburst de- 
pend on the mass of the planet for a fixed a and accretion 
rate Mi n - All properties are derived from the simulations 
apart from i?n m ,LC and ti ast ,LC, which are our semi-analytic 
predictions (see below). Fig.|H]shows the corresponding light 
curves. 

Qualitatively, we expect 7? t ri g to increase with planet 



mass, as the disc has to bank up to higher surface densities 
to push the planet inwards. In p r incipl e, we could use the 
analytic solutions of llvanov et alJ Ill999 j ) to predict fltrfg (i n 
a manner similar to that attempted bv lClarke fc Sverll 996h 
but in practice we found t hat the simulation only converges 
to the llvanov et alJ ( 1999) solution when the planet is about 
a factor 10 inside Rb, the radius at which banking up of 
the disc first becomes significant. In practice, the instability 
is triggered when the planet is a factor of 2 or so within 
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Table 1. Basic outburst features with varying planet mass. See text for the definition of the various quantities. 

M s /Mj upiteT Mout ipeak/^0 Rtrig/Re Rlim/Re -Rl imiLC /-R0 trcc/yr *last/yr *la S t,Lc/y r 
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Figure 8. Light curves obtained by varying the planet mass. 
Dashed line: 15Afj up it C ri solid line: 10Afj up it cr (this is the 
"standard case"), long dashed line: 5Mj up it C r, dot-dashed 
line: 2Mj upitcr . 



Rb for satel l ites in the planetary mass range and so the 
llvanov et alj IJ1999T) solution is not applicable. 

As noted already, the positive dependence of M ou t and, 
consequently, of L pca k on Rtrig (and hence indirectly on the 
planet mass) can be understood in terms of the larger value 
of fi on the upper branch at Ea when -Rtrig is large. Equation 
1171 however overestimates M ou t, since the upward transi- 
tion when the instability is triggered is not vertical in the 
S-/i plane (see also lBL94T) ; in particular the dependence of 
M ou t on -Rtrig is considerably milder than that predicted by 
equation 117> . 

The outermost location in the disc that is reached by 
the instability, Rii m , is an important quantity since it de- 
termines both the outburst duration and the region of the 
disc co ntributing to the enhanced SED during the outburst 
( Ken von et alJ [1988; Ke nvon fc Ha rtmann 1991). We are 
able to predict Rum from the following simple argument. 
The spike in surface density associated with the ionisation 
front (see, for example, Fig. [I] panel c) must evidently con- 
tain enough mass, when the f r ont is at radius R, to attain 
surface density Ea . iLin et al.1 (1985) have argued that the 
front cannot be narrower than H, since otherwise the ra- 
dial pressure gradient would reverse the gradient of specific 
angular momentum in the disc, thus rendering it Rayleigh 
unstable. Thus, a lower limit to the mass in the ionisation 



front at R is AM(fi) = 2ttRT, a (R)H. During the rise of the 
outburst, mass diffuses out to the ionisation front from pre- 
viously destabilised regions at smaller radii. We thus derive 
an estimate of the outermost propagation radius of the in- 
stability, -Rii m ,LC> as being the radius at which AM (Rii m ,Lc) 
is equal to the mass interior to -Rn m ,LC when the outburst 
is first trigge r ed. S ince, as noted above, we cannot use the 
llvanov et alJ J1999T) solution to predict the surface density 
profile in the disc when the outburst is first triggered, we 
instead obtain the interior mass from the simulation. The 
resulting estimator Rii mj LC is in excellent agreement with 
that obtained in the simulation. We also find that the vis- 
cous timescale of the disc on the upper branch at Rii m ,LC 
provides a good estimator of outburst duration (ti as t,Lc), if 
outburst duration is defined as being the period over which 
the luminosity is enhanced by more than a fa ctor 2 compared 
with the pre-outburst level. [We note that BL94 proposed 
a different estimator of Ru m , based on the minimum radius 
at which the disc can remain on the lower stable branch at 
an ac cretion rate of Mm. However, as noted by iBell et all 
Il995l i?n m can exceed this value in the case of a triggered 
outburst.] 

Finally, we note that the recurrence time, t rcc , also in- 
creases with planet mass, since larger planets, for which Ru m 
is larger, are associated with outbursts that clear out the 
disc over a larger radial range. We found, however, that the 
recurrence timescale is not exactly linearly proportional to 
the mass accreted onto the star during the outburst, being 
slightly longer than expected for the less massive planets. 
We attribute this effect to the fact that in order to replenish 
the disc through M- m during quiescence, the matter has to 
diffuse to smaller radii in the case of less massive planets 
(since -Rtrig is smaller), hence slowing down the process. 



4-2.2 Accretion rate 

When the accretion rate Mi n , is lowered by an order of mag- 
nitude to 3 10 _7 MQ/yr the planet spirals in more slowly 
(see Fig. [51 left panel). This longer timescale means that the 
angular momentum that the planet must lose in order to spi- 
ral in to a given radius can be removed at a lower rate, and 
hence involves a less steep banking up of the surface density 
in the disc. Consequently, -Rtrig is attained at smaller radius 
(5R©, compared to 10R© in the reference case. The smaller 
value of Rtrig results in a lower amplitude, shorter duration 
outburst, as expected (Fig. |U] right panel). 



4.2.3 Value of a 

The value of a on the high branch influences the outburst in 
a manner similar to the influence of this parameter on un- 
triggered outbursts. Thus when cVhigh is reduced by a factor 
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Figure 9. Left: Orbital evolution of the planet for two different choices of M\ n = 3 10 Mq/jt (solid line) and 3 10 Mm/yr (dotted 
line). Right: Corresponding light curves. 
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Figure 10. Light curves for two different choices of cthigh'- (solid 
line) «hi g h = 5 10 -4 , (dashed line) cthigh = 10 — 3 . 



two, the rise and decay timescales are somewhat larger (Fig. 
11011 . On the other hand, Eg is increased (Equation[7J so the 
disc has to accrue a smaller increase in surface density to 
trigger the next outburst and consequently the recurrence 
time is decreased with respect to the reference case from 
w 3800 years to « 2000 years. 



5 DISCUSSION 

Thermal instability models for FU Orion i s obje cts have been 
investigated previously bv lClarke et al.l il990f) (who inves- 
tigated outburst artificially triggered in an otherwise stable 
disc), by IBL94 (who considered "self-regulated" outbursts, 
produced in discs fed at a high enou gh rate to be the rmally 
unstable in the inner regions) and bv lBell et alj <I1995T) (who 
investigated artificially triggered outbursts in discs fed at a 
high rate). 

All these previous investigations have used a more re- 
fined energy equation with respect to ours, since they have 
retained both the advective and the radial radiative fluxes 
terms, that we neglect here. The imp ortance of these non - 
local effects has been noted initially bv lClarke et alj |1989), 
in determining th e outer propagation radius of the instabil- 
ity. Subseauentlv. lClarke et ah! il990Tl have noted that, in an 
outside-in outburst, as the instability propagates inwards, it 
needs to heat up the inner disc (which is initially much colder 
than the temperature of partial ionisation of hydrogen) be- 
fore making it jump to the upper branch. Therefore, initially 
most of the energy gained through the enhanced accretion, 
rather than being radiated away, is used to heat up the inner 
disc, resulting in a sharper transition to the high luminosity 
phase (on the thermal timescale of the inner disc) , once the 
front reaches the star. In our simulations, we are not able to 
follow this behaviour, since we do not include any advective 
term in our energy equation. We expect that the inclusion of 
these terms will result in even shorter rise times compared 
to the iriso ~ 2 yrs that we find here. 

The outbursts have also interesting consequences for 
the planet migr ation. We have in fa ct confirmed the pre- 
diction made bv lClarke fc Sverl (119961) that the migration of 
the planet is going to be significantly slowed down after the 
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planet has reached the radius at which the instability is first 
triggered. 

If, on the one hand, we have demonstrated how the trig- 
gering effect due to an embedded planet is able to reduce 
the rise timescale, which in untriggered models would be 
too large compared to observations, the duration of the out- 
burst still remains an issue for all thermal instability models. 
As we discussed above, the duration of the outburst is re- 
lated to the outermost propagation radius of the instability, 
which, in thermal instability models, is bound to be rela- 
tively small (see discussion in Sec. l4.2H . Even in our longest 
lived outburst, £i ast was not larger than ~ 60 years, whereas, 
for example, the duration of the outburst in FU Ori itself 
appears to be larger than 70 years. In order to reproduce 
these long durations we needed to stretch our parameters 
to rather unrealistic values, with very high Mj n (larger than 
5 lO- 5 A-/ /yr) and low a Mgh (« 5 10~ 4 ). 

Our simplified model can be refined in several ways: 
of course, to obtain more realistic light curves, we should 
add the neglected terms in the energy equation, in partic- 
ular the advective terms, that might influence the shape of 
the light curve, as discussed above. A second limitation of 
our model is that we have simply assumed that the torque 
between the planet and the disc is given by equation @. 
As already discussed in Sectional while this assumption is 
reasonable in the quiescent phase of the outburst, when the 
planet opens up a gap in the disc, its use during the outburst 
phase is questionable, since during the outburst the disc is 
going to fill the gap and most of the torque exerted on the 
planet might come from a circumplanetary disc. In addi- 
tion, corotation torques mig ht also play a role in this case 
iJMasset fc Papaloizoul 120031) . Most 2D and 3D numerical 
simulations of this process (iLu bow et all H 999: Bate et alJ 
12003 : iMasset fc Papaloizoul I2003T) have usually considered 
the case of a planet embedded in a much colder, T Tauri-like, 
disc and are therefore not appropriate to test the behaviour 
of the planet-disc system in the case we are interested in 
here. Therefore, it would be important to perform detailed 
simulations of planet-disc interaction in FU Orionis objects. 

A full study of these effects is highly desirable because 
the early evolution of the planet is going to be strongly in- 
fluenced by the outbursts, both because the detailed planet- 
disc interaction is going to affect the planet migration and 
the mass accretion rate onto the planet itself (note that in 
the present work the planet mass is fixed) and because we 
expect that the structure of the planet will be strongly af- 
fected by the very hot "thermal bath" in which the planet 
is embedded during the outburst. 

From a more general point of view, we note that the 
process described here implies that planets are present in 
the early stages of star formation, during which FU Orio- 
nis outbursts are thought to occur. On the other hand, core 
accretion models of planet formation (Pol lack et alJ 11996) 
would imply very long timescales for the formation of gas 
giant planets (this actually is one problem of the core ac- 
cretion scenario itself). Assessing the likeliness of forming 
a planet early in the star formation process is beyond the 
scopes of this paper and only a full theory of planet for- 
mation will give some clues. However, we point out that (i) 
the planets considered here are probably going to be swal- 
lowed by the star after a few outbursts, and could therefore 
represent an early generation of planets with only a minor 



fraction of them surviving the pre-main-sequence evolution 
of the star; (it) we are here considering high mass planets, 
which are not easy to form via core accretion and might have 
formed through some other proces s, such as g ravitational in- 
stabilities in the protostellar disc (Boss 2000). Finally, (Hi), 
we note that the youth of FU Orionis objects has been gen- 
erally inferre d from the high input mass accretion rates re- 
quired in the|BL94] model to trigger the thermal instability. 
Here we have shown that a planet-induced outburst may 
also occur with a lower Mj n , which would allow FU Orio- 
nis outbursts to occur at a later stage in the protostellar 
evolution. 



6 CONCLUSIONS 

In this paper we have shown how the interaction between a 
protostellar disc and an embedded massive planet is going 
to naturally produce outbursts in the disc, due to the bank- 
ing up of material upstream of the planet position, which 
triggers a thermal instability in the disc. The location at 
which the instability is triggered depends on the mass of 
the embedded planet, and is larger for more massive plan- 
ets. For a lOAf jupitor planet, we obtain i?tri g ~ 107?©. In this 
way, the instability propagates outside-in, therefore produc- 
ing rapid-rise outbursts (with i r isc « 2 yrs), as observed in 
some FU Orionis objects (such as FU Ori itself and V1057 
Cyg). Longer rise times, as observed in the case of V1515 
Cyg, can be obtained either in non-triggered outbursts, of 
the kind discussed by |BL94 . or in systems with a low-mass 
embedded planet (with M s < 2 Mj UP it er). We have also revis- 
ited the arguments often used <lBL94T) to estimate the outer 
propagation radius of the instability, Ru m , since these pre- 
vious arguments are not appropriate in the case of triggered 
outbursts, such as those described in this paper. 

The occurrence of FU Orionis outbursts can have im- 
portant effects on the embedded planet, (i) from a dynami- 
cal point of view, since we have demonstrated here that the 
outburst mechanism is going to significantly slow down the 
inward migration of the planet, which is going to linger near 
the position where the first outburst is triggered (see Fig. 
EJ, and (it) concerning the internal structure of the planet, 
which as a consequence of the outburst will be suddenly em- 
bedded in a very hot environment. We have not discussed 
this last topic in the present work and we leave it to further 
investigations. 
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